##########################################################################################################################################
###### This is the file containing the regression results for the tail risk analysis
###############################
# It creates the final results for Table 8 (and part of Table 1)
# Please note that in some instances, the paper contains some formatting modifications
# such as a "log()" that is not automatically added here or that we present the negative
# of tsim, while most output here is generated for the positive

# This is just the final step. There are other code files that are used for data preparation
# This processed data are saved and loaded in this code

# set to your working directory:
setwd("C:\\Users\\Admin\\OneDrive\\Desktop\\Market Power and Systematic Risk - Data and Code")
getwd()


######
# this is a file containing the tail risk estimates for individual companies
load("BT11Q30.RData")
#####
# this is a file containing information on the companies in the BT11Q30 file
SECID = read.table("Info_Company.txt")
#Write The Info_Company differently
write.table(SECID,file="SECID.txt",
            row.names=F)

for(i in 1:length(finalMatrix)){
  Temp = finalMatrix[[i]]
  if(is.null(nrow(Temp))){
    
  }else{
    finalMatrix[[i]]$SECID = rep(SECID[i,1],nrow(Temp))
  }
  print(i)
}

for(i in length(finalMatrix):1){
  if(is.null(nrow(finalMatrix[[i]]))|!is.null(finalMatrix[[i]]$message)){
    finalMatrix[[i]] = NULL
  }else if(nrow(finalMatrix[[i]])<=1){
    finalMatrix[[i]] = NULL
  }
}

library(data.table)
TailRisks <- rbindlist(finalMatrix)
TailRisks$YearMonth = paste0(format(TailRisks$Date,"%Y"),format(TailRisks$Date,"%m"))

Temp = setDT(TailRisks)[order(Date), lapply(.SD,mean,na.rm=T), by = c("YearMonth","SECID"),.SDcols=c("BT11Q_LT","BT11Q_RT","TLM")]
Temp_Two = setDT(TailRisks)[order(Date), .SD[which.max(Date)], by = c("YearMonth","SECID"),.SDcols=c("Date")]
TailRisks = merge(Temp,Temp_Two,by=c("YearMonth","SECID"))

######
# this is a file containing data from OptionMetrics to match the secid identifiers
# to others. The file contains the following elements from OptionMetrics (in that order):
# - secid	
# - ticker
# - cusip 
# - sic
# - effect_date
Matching = read.table("SECMatching.txt",sep="\t",header=T)
#####

Matching$effect_date = as.Date(as.character(Matching$effect_date),format="%Y%m%d")

IDS <- unique(TailRisks$SECID)
TailRisks$cusip = rep(NA,nrow(TailRisks))
options(warn=2)
setDF(TailRisks)
for(i in 1:length(IDS)){
  TempSemi = TailRisks[TailRisks$SECID==IDS[i],]
  SubMatch = Matching[Matching$secid == IDS[i],]
  if(nrow(SubMatch)>=1){
    for(j in 1:nrow(SubMatch)){
      TempSemi$cusip[TempSemi$Date>SubMatch$effect_date[j]] = rep(SubMatch$cusip[j],length(TempSemi$cusip[TempSemi$Date>SubMatch$effect_date[j]]))
    }
    TailRisks[TailRisks$SECID==IDS[i],] = TempSemi
  }
  print(i)
}
options(warn=1)

load("Control_Panel.RData")
setDT(TailRisks)
TailRisks <- na.omit(TailRisks)
table(TailRisks[,length(unique(cusip)),by="SECID"]$V1)
Double_Cusip= TailRisks[,length(unique(cusip)),by="SECID"]$SECID[which(TailRisks[,length(unique(cusip)),by="SECID"]$V1>1)]
Test_vektor = c()
for(i in 1:length(Double_Cusip)){
  Test_vektor[i] = all(unique(TailRisks[TailRisks$SECID == Double_Cusip[i],]$cusip)%in%Panel_New$CUSIP)
}

Test_vektor = c()
for(i in 1:length(Double_Cusip)){
  Test_vektor[i] = any(unique(TailRisks[TailRisks$SECID == Double_Cusip[i],]$cusip)%in%Panel_New$CUSIP)
}

Single_Cusip= TailRisks[,length(unique(cusip)),by="SECID"]$SECID[which(TailRisks[,length(unique(cusip)),by="SECID"]$V1==1)]
Test_vektor = c()
for(i in 1:length(Single_Cusip)){
  Test_vektor[i] = all(unique(TailRisks[TailRisks$SECID == Single_Cusip[i],]$cusip)%in%Panel_New$CUSIP)
}

Panel_New<- merge(Panel_New,TailRisks,by.x=c("CUSIP","YearMonth"),by.y=c("cusip","YearMonth"))


Temp = Panel_New[,c("Time","YearMonth","year","NAICS","PERMNO","BT11Q_LT","BT11Q_RT","TLM",
                    "HHI","prodmktfluid","tnic3hhi","tnic3tsimm","Age","Momentum","Illiquidity","MktCap","Mkt_Book","Default_Spread","Financial_Leverage", "Operating_Leverage","Firm_size","InvestmentRates","ROE","Logat","Dividend_Dummy","q","I_Vol")]
for(i in 1:dim(Temp)[2]){
  Temp[,i][Temp[,i] == Inf] = NA
  Temp[,i][Temp[,i] == -Inf] = NA
}
Temp = na.omit(Temp)

library(lfe)

source("Table_Functions.R")

Temp$HHI=scale(log(Temp$HHI))
a<-poly(Temp$HHI,2)
Temp$HHI2=a[,2]/sd(a[,2])
Temp$prodmktfluid=scale(Temp$prodmktfluid)
a<-poly(Temp$prodmktfluid,2)
Temp$prodmktfluid2=a[,2]/sd(a[,2])
Temp$tnic3hhi=scale(Temp$tnic3hhi)
a<-poly(Temp$tnic3hhi,2)
Temp$tnic3hhi2=a[,2]/sd(a[,2])
Temp$tnic3tsimm=scale(Temp$tnic3tsimm)
a<-poly(Temp$tnic3tsimm,2)
Temp$tnic3tsimm2=a[,2]/sd(a[,2])
Temp$Age=scale(log(Temp$Age))
Temp$Momentum=scale(Temp$Momentum)
Temp$Illiquidity=scale(Temp$Illiquidity)
Temp$MktCap=scale(log(Temp$MktCap))
Temp$Mkt_Book=scale(log(Temp$Mkt_Book))
Temp$Financial_Leverage=scale(Temp$Financial_Leverage)
Temp$Operating_Leverage=scale(Temp$Operating_Leverage)
Temp$Default_Spread=scale(Temp$Default_Spread)
Temp$Firm_size=scale(Temp$Firm_size)
Temp$InvestmentRates=scale(Temp$InvestmentRates)
Temp$ROE=scale(Temp$ROE)
Temp$Logat=scale(Temp$Logat)
Temp$q=scale(Temp$q)
Temp$I_Vol=scale(Temp$I_Vol)

LT_prod = felm(BT11Q_LT ~ tnic3tsimm + tnic3tsimm2 + HHI + HHI2 + ((Age)) + (Logat) + (Default_Spread) + (Dividend_Dummy) + (Financial_Leverage) + 
                   ((MktCap)) + (Illiquidity) + (InvestmentRates) + (I_Vol) + ((Mkt_Book)) + (Momentum) +
                   (Operating_Leverage) + (q) + (ROE) | PERMNO + year | 0 | NAICS + year, data = Temp, cmethod="reghdfe",psdef=FALSE)

library(fixest)
model_feols <- feols(data = Temp, BT11Q_LT ~ tnic3tsimm + I(tnic3tsimm^2) + HHI + HHI2 + ((Age)) + (Logat) + (Default_Spread) + (Dividend_Dummy) + (Financial_Leverage) + 
                       ((MktCap)) + (Illiquidity) + (InvestmentRates) + (I_Vol) + ((Mkt_Book)) + (Momentum) +
                       (Operating_Leverage) + (q) + (ROE) | PERMNO + year)
newdata=Temp[1:5,]
newdata[2,]=newdata[1,]
newdata[3,]=newdata[1,]
newdata[4,]=newdata[1,]
newdata[5,]=newdata[1,]
newdata$tnic3tsimm[1]<--2
newdata$tnic3tsimm[2]<--1
newdata$tnic3tsimm[3]<-0
newdata$tnic3tsimm[4]<-1
newdata$tnic3tsimm[5]<-2
predict(model_feols, newdata)




RT_prod = felm(BT11Q_RT ~ tnic3tsimm + tnic3tsimm2 + HHI + HHI2 + ((Age)) + (Logat) + (Default_Spread) + (Dividend_Dummy) + (Financial_Leverage) + 
                  ((MktCap)) + (Illiquidity) + (InvestmentRates) + (I_Vol) + ((Mkt_Book)) + (Momentum) +
                  (Operating_Leverage) + (q) + (ROE) | PERMNO + year | 0 | NAICS + year, data = Temp, cmethod="reghdfe",psdef=FALSE)

Partial_Paper = data.frame(Coeff_Names = rep(NA,2*length(coef(LT_prod))+2),LT = rep(NA,2*length(coef(LT_prod))+2), 
                           RT = rep(NA,2*length(coef(LT_prod))+2))

Names = c("tsim","tsim^2","log(HHI)","log(HHI)^2","log(Age)","log(AT)","Default~spread","Dividend","Financial~leverage","log(Firm~size)","Illiquidity","Investment~rate","iVol",
          "log(Mkt/Book)","Momentum","Operating~leverage","q","ROE")


Partial_Paper$Coeff_Names[seq(from=1,to=(nrow(Partial_Paper)-2),by=2)] = Names
Partial_Paper$Coeff_Names[nrow(Partial_Paper)-1] = "R^2"
Partial_Paper$Coeff_Names[nrow(Partial_Paper)] = "Fixed-Effects"

Partial_Paper$LT[seq(from=1,to=(nrow(Partial_Paper)-2),by=2)] = paste0(as.character(Lag_zero(round(summary(LT_prod)$coefficients[,1],digits=3),digits=3)),star_function(summary(LT_prod)$coefficients[,4]))
Partial_Paper$LT[seq(from=2,to=(nrow(Partial_Paper)-2),by=2)] = paste0("(",Lag_zero(round(summary(LT_prod)$coefficients[,3],digits=3),digits=3),")")
Partial_Paper$LT[nrow(Partial_Paper)-1] = as.character(Lag_zero(round(summary(LT_prod)$r2adj,digits=4) * 100,digits=2))
Partial_Paper$LT[nrow(Partial_Paper)] = "Yes"

Partial_Paper$RT[seq(from=1,to=(nrow(Partial_Paper)-2),by=2)] = paste0(as.character(Lag_zero(round(summary(RT_prod)$coefficients[,1],digits=3),digits=3)),star_function(summary(RT_prod)$coefficients[,4]))
Partial_Paper$RT[seq(from=2,to=(nrow(Partial_Paper)-2),by=2)] = paste0("(",Lag_zero(round(summary(RT_prod)$coefficients[,3],digits=3),digits=3),")")
Partial_Paper$RT[nrow(Partial_Paper)-1] = as.character(Lag_zero(round(summary(RT_prod)$r2adj,digits=4) * 100,digits=2))
Partial_Paper$RT[nrow(Partial_Paper)] = "Yes"

library(xtable)
print(xtable(Partial_Paper),include.rownames=F,sanitize.text.function=function(x){x})


#######################################
# Table 1: Summary Statistics (Lines 8-9)


setDT(Temp)
mean(Temp[,mean(BT11Q_LT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,mean(BT11Q_RT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,sd(BT11Q_LT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,sd(BT11Q_RT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,median(BT11Q_LT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,median(BT11Q_RT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,min(BT11Q_LT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,min(BT11Q_RT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,max(BT11Q_LT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,max(BT11Q_RT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,skewness(BT11Q_LT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,skewness(BT11Q_RT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,kurtosis(BT11Q_LT,na.rm=T),by=YearMonth]$V1,na.rm=T)
mean(Temp[,kurtosis(BT11Q_RT,na.rm=T),by=YearMonth]$V1,na.rm=T)
setDF(Temp)

############################

